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Preface 


Earlier  in  this  study,  the  emphasis  was  on  under¬ 
standing  the  occurrence  of  oscillatory  instabilities 
in  finite  difference  approximations  for  two-dimensional 
transient  heat  conduction (diffusion) .  The  Peaceman- 
Rachford  alternating  direction  implicit  (ADI)  finite 
difference  method  (FDMTH)  was  of  special  interest. 

As  the  computer  programs  for  the  various  FDMTH s  were 
debugged  and  data  began  to  accumulate,  the  emphasis 
shifted  to  the  accuracy  of  the  FDMTHs .  The  Crank - 
Nicholson  implicit  FDMTH  proved  to  be  the  most  accurate 
of  the  methods  considered  and  the  Peaceman-Rachford  ADI 
FDMTH  the  least  accurate. 

Dr.  Bernard  Kaplan's  guidance  and  encouragement 
throughout  this  study  were  always  timely  and  effective. 
Special  thanks  are  due  Dr.  W.  Kessler  of  the  Air  Force 
Materials  Laboratory  for  sponsoring  this  research  project. 
I  am  especially  indebted  to  my  wife,  Madalene,  and  our 
three  children  who  have  survived  my  master's  thesis. 


T.  Sidney  Chivers,  Jr. 
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Abstract 


\ 

~  The  two-dimensional  transient  heat  conduction 
(diffusion)  equation  was  solved  using  the  fully  explicit, 
fully  implicit,  Crank-Nicholson  implicit,  and  Peaceman- 
Rachford  alternating  direction  implicit  (ADI)  finite 
difference  methods  (FDMTHs) .  The  general  stability 
condition  for  the  same  FDMTHs  was  derived  by  the  matrix, 
coefficient,  and  a  probabilistic  method.  The  matrix, 
coefficient,  and  probabilistic  methods  were  found  to  be 
equivalent  in  that  each  lead  to  the  same  general  stability 
condition.  Oscillatory  behavior  of  the  fully  explicit 
FDMTH  was  as  predicted  by  the  general  stability  condition. 
Though  the  Crank-Nicholson  implicit  and  the  Peaceman- 
Rachford  ADI  FDMTHs  were  expected  to  be  unconditionally 
stable,  unstable  oscillations  were  observed  for  large 
sizes  of  time  step.  For  large  numbers  of  time  steps 
and  sizes  of  time  steps  for  which  all  FDMTHs  considered 
are  stable,  the  Crank-Nicholson  implicit  FDMTH  is  the 


more  accurate 


I .  Introduction 


Background . 


^3  Finite  difference  methods  are  useful  in  obtaining 
solutions  for  engineering  problems  involving  partial 
differential  equations  that  cannot  be  solved  in  closed 
form.  Because  finite  difference  methods  approximate 
the  true  solution,  their  competent  use  requires  an 
understanding  of  discretization  errors  and  the  stability 


condition.  The  discretization  error  is  the  combined 
effect  of  round-off  and  truncation  due  to  the  "limitation 
on  the  number  of  significant  figures  carried  by  a 
computer"  (2s 20)  and  the  truncation  of  higher  order 
Taylor  series  terms  in  developing  the  finite  difference 
approximations  of  partial  differentials  (2:20), 
respectively.  The  stability  condition  defines  parameter 
regions  in  which  the  finite  difference  method  remains 
stable  for  large  numbers  of  time  steps.  The  coefficient 
method  (15,16,12:283),  the  matrix  method  (20:60-68),  and 
the  Fourier  method  (13)  are  the  more  common  methods  of 
deriving  the  stability  condition.  A  probabilistic 
method  of  deriving  the  stability  condition  is  suggested 
by  the  work  of  Kaplan  (10) . 


Problem  Statement. 


The.  primary  objectives  of  this  study  were  to 
understand  instability  in  finite  difference  methods 
and  be  able  to  predict  when  oscillatory  behavior  would 
occur.  Secondary  objectives  were  to  compare  the  various 
methods  of  deriving  the  stability  condition,  and  compare 
finite  difference  methods  on  the  basis  of  discretization 
error  and  stability. 

Scope ♦ 

This  study  was  limited  to  two-dimensional  transient 
heat  conduction (diffusion)  in  a  rectangular  region. 

The  finite  difference  methods  considered  were  the 
fully  explicit,  fully  implicit.  Crank -Nicholson  implicit 
and  Peaceman-Rachford  alternating  direction  implicit. 

The  matrix,  coefficient,  and  probabilistic  methods  of 
deriving  the  stability  condition  were  considered. 


Computer  programs  were  developed  to  solve  the 
two-dimensional  transient  heat  conduction (diffusion) 
problem  by  either  the  fully  explicit,  fully  implicit, 
Crank-Nicholson  implicit,  or  Peaceman-Rachford 
alternating  direction  implicit  finite  difference  method. 
The  thermodynamic  and  mathematical  aspects  of  the 
stability  of  finite  difference  methods  were  researched. 
The  stability  condition  for  the  general  two-dimensional 
finite  difference  approximation  of  transient  heat 
conduction (diffusion)  was  derived  by  the  matrix, 
coefficient,  and  probabilistic  methods.  Computer 
programs  developed  were  run  for  selected  time  increments 
and  nodal  array  sizes  to  develop  data  for  comparison 
of  finite  difference  methods,  and  to  assess  the 
validity  of  the  stability  condition  derived  for  the 
general  finite  difference  method. 


Sequence  of  Presentation. 

Finite  difference  methods  are  discussed  in  the 
next  section.  The  matrix  methods  used  are  described 
in  Section  III.  The  theory  of  stability  analysis  is 
detailed  in  Section  IV.  Computer  methods  and  other 
procedures  used  are  described  in  Section  V.  Results 
are  summarized  in  Section  VI  with  graphic  results  appended. 

Acronyms . 

Two  acronyms  will  be  used  throughout  the  remainder 
of  this  thesis.  FDMTH(s)  will  represent  finite  difference 
method (s),  and  ADI  will  represent  alternating  direction 
implicit.  These  acronyms  are  necessary  for  brevity. 

Using  these  acronyms,  the  Peaceman-Rachford  alternating 
direction  implicit  finite  difference  method  becomes  the 


Peaceman-Rachford  ADI  FDMTH. 


Transient  Heat  Conduction (Diffusion) 


by  Finite  Differences 


The  two-dimensional  transient  heat  conduction 
(diffusion)  problem  considered  in  this  study  is  described 
by  the  following  initial-boundary  value  problem. 


3T 

3  2T 

32T 

=  a  I  - 

+  -  J 

(la) 

3t 

3x2 

3y2 

t) 

=  T  (x ,  0 ,  t) 
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where  T  is  for  temperature,  a  is  the  thermal  diffusivity, 
x  and  y  are  spatial  dimensions,  and  t  is  for  time.  The 
analytic  solution,  derived  at  Appendix  A,  is 

T(x,y,t)  =  exp[  -a(m2TT2  +  n2ir2)t  ] 

sin(mirx)  sin(niry)  (Id) 

In  the  following,  the  two-dimensional  finite  difference 
approximation  of  transient  heat  conduction (diffusion) 
is  developed  from  two  basic  difference  equations,  the 
first  forward  difference  and  the  second  central  difference 
Error  analysis  is  discussed  in  the  final  sub-section. 


The  General  Finite  Difference  Approximation 


Approximating  the  time  derivative  in  equation  (1) 
by  a  first  forward  difference  and  the  two  spatial 
derivatives  by  second  central  differences,  the  finite 
difference  approximation  of  two-dimensional  transient 
heat  conduction (diffusion)  becomes 


T(i,j,k+1)  -  T(i,j,k) 


o [  T(i+1,  j)  -  2T(i,j)  +  T(i-l,j)  ] 


a  l  T (i, j+1)  -  2T(i,j)  +  T (i, j-1)  ] 


where  the  subscripts  i,  j,  and  k  are  used  to  identify 
the  node  location  and  time  step.  Figure  1  depicts  a 
spatial  domain  discretized  by  imposing  a  5  by  5  array 
of  nodes  over  the  domain  so  that,  for  equation  (2) , 
i  =  0,1,2,  ...  5-1  and  j  =  0,1,2,  ...  5-1.  For  a 
general  m  by  n  nodal  array,  i  =  0,1,2,  ...  m-1  and 


Figure  2.  Relative  Labelling 
of  Nodes 


Equation  (2)  now  simplifies  to 
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The  temperatures  on  the  right  side  of  equation  (2) 
represent  a  mean  temperature  between  time  steps  (15:28) 
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where  f  is  some  weighting  factor. 


From  equations  (8)  and  (9)  ,  the  general  two-dimensional 
FDMTH  for  transient  heat  conduction (diffusion)  becomes 
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Figure  3  indicates  the  relative  locations  of  the 

temperatures  in  equation  (10) .  Values  of  f  can  be 

e 

varied  to  reduce  equation  (10)  to  a  specific  finite 

difference  method.  Table  1  lists  values  of  f  for 

e 

the  fully  explicit,  fully  implicit,  Crank-Nicholson 
implicit,  and  Peaceman-Rachford  ADI  FDMTH s . 


(10) 


Table  1 


Weighting  Factors  for  the  Generalized 
Finite  Difference  Approximation  of 
Transient  Heat  Conduction  (Diffusion) 


FDMTH 

fE 

£w 

fp 

EW 

£n 

£s 

fp 

*NS 

Fully  Explicit 

0 

0 

0 

0 

0 

0 

Fully  Implicit 

1 

1 

1 

1 

1 

1 

Crank -Nicholson 

Implicit 

>5 

% 

H 

% 

*5 

Peaceman-Rachford  ADI 

for  Odd  Time  Steps 

0 

0 

0 

1 

1 

1 

for  Even  Time  Steps 

1 

1 

l 

0 

0 

0 
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Fully  Explicit  FDMTH. 


The  two-dimensional  fully  explicit  FDMTH  is 
explicit  in  both  spatial  dimensions.  This  FDMTH  assumes 
nodal  temperatures  for  time  t  prevail  until  time  t  +  At 
(15:56).  From  equation  (10)  and  Table  1,  the  equation 
for  the  fully  explicit  FDMTH  for  transient  heat 
conduction (diffusion)  is 
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Figure  4  is  the  mnemonic  for  the  fully  explicit  FDMTH. 
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Fully  Implicit  FDMTH . 

The  two-dimensional  fully  implicit  FDMTH  is 
implicit  in  both  spatial  dimensions.  This  FDMTH  assumes 
nodal  temperatures  at  time  t  change  immediately  to 
temperatures  for  time  t  +  At  which  prevail  throughout 
the  time  step.  From  equation  (10)  and  Table  1,  the 
equation  for  the  fully  implicit  FDMTH  for  transient 
heat  conduction (diffusion)  is 
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Figure  5  is  t^e  mnemonic  for  the  fully  implicit  FDMTH. 
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i'igure  5.  Mnemonic  for  the  Fully  Implicit 
FDMTH. 
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Crank-Nicholson  Implicit  FDMTH. 


The  two-dimensional  Crank-Nicholson  FDMTH  assumes 
nodal  temperatures  change  linearly  from  their  value  at 
time  t  to  their  value  at  time  t  +  At.  From  equation 
(10)  and  Table  1,  the  equation  for  the  Crank-Nicholson 
implicit  FDMTH  for  transient  heat  conduction (diffusion) 


1  +  - j  +  - 1 

(Ax)  2  (Ay) 2 


|mk+l  . 

2  'e  W  ' 
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UN  XS  1 
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Figure  6  is  the  mnemonic  for  the  Crank-Nicholson 
implicit  FDMTH. 
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Peaceman-Rachford  ADI  FDMTH 


The  Peaceman-Rachford  ADI  FDMTH  is  a  two  step 
method.  In  the  first  step,  nodal  temperature  changes 
are  implicit  with  respect  to  one  spatial  dimension  and 
explicit  with  respect  to  the  other.  In  the  second 
time  step,  the  explicit /implicit  roles  of  the  two 
spatial  dimensions  are  reversed.  From  equation  (10) 
and  Table  1,  the  equations  for  the  Peaceman-Rachford 
ADI  FDMTH  for  transient  heat  conduction (diffusion)  are 
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Figures  7  and  8  are  the  mnemonics  for  the  Peaceman- 


Rachford  ADI  FDMTH 


r  r  r  r 


Error  Analysis. 


The  discretization  error  is  the  sum  of  the  round-off 
error  and  the  truncation  error.  The  round-off  error  is 
dependent  on  the  largest  number  of  digits  that  can  be 
represented  in  a  computer's  memory  and  the  number  of 
computations  needed  to  obtain  a  solution.  The  truncation 
error  is  due  to  the  truncation  of  higher  order  terms  of 
the  Taylor  series  representation  of  differentials  in 
approximating  the  differentials  by  finite  differences 
(2:20).  Table  2  lists  the  order  of  truncation  error 
expected  with  the  fully  explicit,  fully  implicit, 
Crank-Nicholson  implicit,  and  Peaceman-Rachford  ADI 
FDMTHs .  In  this  study,  the  discretization  error,  ERR, 
is  computed  as  the  difference  between  the  FDMTH  and 
analytic  solution  temperatures,  or 


k+1 

where  T~  is  the  FDMTH  temperature  of  node  e  at  time 
k+1 

step  k+1  and  TRv  is  the  analytic  solution  temperature 
of  node  e  at  time  step  k+1. 

Numerical  stability,  another  aspect  of  error  analysis 


of  FDMTHs,  is  discussed  in  Section  IV. 


Table  2 


Order  of  Truncation  Error  of 
Selected  Finite  Difference  Methods 


FDMTH 

Fully  Explicit 

Fully  Implicit 

Crank -Nicholson 

Implicit 

Peaceman-Rachford  ADI 

for  Odd  Time  Steps 


Order  of  the 
Truncation  Error 

At  +  (Ax)2  +  (Ay)2 

At  +  (Ax)2  +  (Ay)2 

(At)2  +  (Ax)2  +  (Ay) 

At  +  (Ax)2 
At  +  (Ay)2 


for  Even  Time  Steps 


Ill .  Matrix  Methods 

Matrix  methods  are  used  to  solve  the  simultaneous 
equations  that  result  when  applying  any  FDMTH  that  is 
partially  or  fully  implicit/  such  as  the  fully  implicit 
Crank-Nicholson  implicit,  and  Peaceman-Rachford  ADI 
FDMTHs.  The  matrix  equivalent  of  the  general  FDMTH 
equation,  equation  (10) ,  is 

A  Tk+1  =  Tk  (19 

where  the  coefficient  matrix  A  is  as  shown  in  Figure  9. 
For  example,  the  coefficient  matrix  for  a  5  by  5  nodal 
array  imposed  over  a  square  domain  for  the  initial¬ 
boundary  value  problem  (1) ,  may  be  represented  as 


Implicit  FDMTH 


where  the  elements  not  shown  are  zero  valued  and 


a  =  1  +  2oAt[l-f_  J  /  ( Ax)  2  +  2aAt  [1-fp  ]/(Ay)2 


EW 


NS 


”  al,l  “  a2 , 2  *  **•  =  an,n 


(21) 


b  2  aAt(l-fx)/(Ax)z  ,  f x  =  f E  =  f w 


“  al,2  “  a2,3  *  ***  “  an-l ,n 


“  a2,l  ~  a3 / 2  *  **•  =  anf n-1 


(22) 


c  =  aAt(l-fy)/(Ay)z  ,  fy  =  fN  =  fs 


al,M  a2,M+l  ~  an-M,n 


aM,l  aM+l ,2 


a 


n ,  n-M 


(23) 


These  type  matrices  are  referred  to  as  banded  because 
all  non-zero  elements  are  within  a  band  of  n  diagonals 
centered  on  the  principal  diagonal  of  the  coefficient 
matrix.  The  width  of  the  band  is  dependent  on  the  array 
size  used  and  the  order  in  which  the  interior  nodes  are 
numbered.  For  example,  if  the  interior  nodes  of  a  8  by  4 
nodal  array  are  numbered  from  left  to  right  and  top  to 
bottom,  as  in  Figure  10,  then  the  coefficient  matrix  is  13 
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If  the  interior  nodes  of  the  same  array  are  numbered  top 
to  bottom  and  left  to  right  as  in  figure  11,  then  the 
bandwidth  is  5.  In  general,  the  bandwidth,  BW  can  be 
defined  as 

BW  =  2m  -  3  (24) 

for  an  m  by  n  nodal  array  if  the  interior  nodes  are  first 
numbered  along  the  dimension  corresponding  to  m. 


Figure  11.  Interior  Nodes  of  a 
8  by  4  Nodal  Array 
Numbered  Top  to  Bottom 


If  the  coefficients  in  the  two  outer  diagonals  are  all 
zero,  as  is  the  case  for  the  Peaceman-Rachford  ADI  FDMTH 
then  the  matrix  is  tridiagonal  and  a  solution  can  be 
obtained  by  the  Thomas  method  (3:46-48) .  An  iterative 
process  such  as  the  method  of  Guass-Sidel  (18:40-41) 
will  solve  the  general  form  of  the  coefficient  matrix. 
The  Givens-Householder  transformation  can  be  used  to 
reduce  the  general  form  to  a  tridiagonal  coefficient 
matrix  (14:85-115,22:901-914). 


IV.  Theory  of  Stability  Analysis 


Instability  in  finite  difference  methods  can  be 
considered  caused  by  either  violation  of  the  first  or 
second  laws  of  thermodynamics  (4) ,  or  use  of  an  unstable 
numerical  process.  This  study  concentrates  on  the 
mathematical  stability  of  FDMTHs .  Oscillations  of  the 
same  order  of  magnitude  as  the  true  solution  are  regarded 
as  phenomenon  characterizing  marginal  stability  of  FDMTHs. 
In  the  following,  the  matrix  method  (20:60-68),  coefficient 
method  (15,16,12:281-283),  and  the  probabilistic  method 
of  deriving  the  stability  condition  will  be  described  and 
two-dimensional  forms  presented. 

Matrix  Method. 

A  matrix  equation  for  the  general  FDMTH  for  two- 
dimensional  transient  heat  conduction (diffusion)  is 

A  Tk+1  =  A'  Tk  (28) 

I 

where  the  coefficient  matrices  A  and  A  are  as  described 
in  Appendix  A.  Equation  (28)  can  be  rewritten  as 

[I+(A-I) ]Tk+1  =  [I+(A'-I)]Tk  (29) 

The  stability  condition  is  determined  by  examination  of 
the  eigenvalues  of  [I+(A'-^)J.  The  general  FDMTH  stability 


28 


condition  for  t.fO-dimensional  transient  heat  conduction 
(diffusion) ,  derived  at  Appendix  A,  is 

(1-f  )  (1-f  ) 

k  =  aAt [  - 5-  +  — — 1  <  h  (30) 

(Ax)2  (Ay)2 

where  f x  =  f£  =  fw  and  f  =  fN  =  fg  .  Table  1  depicts 

the  values  of  f  for  the  FDMTHs  studied, 
e 

The  Coefficient  Method. 

Descriptions  of  the  coefficient  method  by  Meyer 

(12:281-283)  and  by  Patankar  (15,16)  differ  only  in 

k+1  k 

symbology.  Each  considers  the  ratio,  X,  of  T  to  T 
as  an  indicator  of  stability.  When  X  is  negative,  the 
FDMTH  is  unstable.  Table  3  gives  the  ranges  of  X 
associated  with  the  four  types  of  stability  of  FDMTHs. 
Figure  12  depicts  the  finite  difference  stability  curves 
for  transient  heat  conduction (diffusion)  in  a  square 
region  with  one  interior  node.  The  stability  condition 
derived  by  the  coefficient  method  is  the  same  as  the 
stability  condition  derived  by  the  matrix  method 
(See  Appendix  B) . 


Table  3 

X  for  the  Four  Possible  Types  of  Behavior  of  FDMTHs 


for  Transient  Heat  Conduction (Diffusion) 
in  a  Rectangular  Region  with  One  Interior  Node 

(20:282) 


k+1 

A>1  Steady,  unbounded  growth.  T  has 

k 

the  same  sign  as  T  and  is  larger 
in  magnitude. 


k+1 

1>A>0  Steady  decay.  T  has  the  same  sign 

]r 

as  T  and  is  smaller  in  magnitude. 

k+1 

0>A>-1  Stable  oscillations.  T  has  the 

Jr 

opposite  sign  as  T  and  is  smaller  in 
magnitude . 

k+1 

\>-l  Unstable  oscillations.  T  has  the 

V 

opposite  sign  as  T  and  is  larger  in 
magnitude . 


Fully  Implicit 


12.  Finite  Difference  Stability  Curves 
for  Transient  Heat  Conduction 
(Diffusion)  in  a  Square  Region  with 
One  Interior  Node 


The  Probabilistic  Method. 


% 


i 


* 
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The  probabilistic  method  is  suggested  by  the  works 
of  Kaplan  (10)  ,  Haji-Sheikh  and  Sparrow  (8) ,  and  Collins 
(5) .  By  considering  the  coefficients  of  the  general 
FDMTH  to  be  probabilistic,  a  stability  condition  can  be 
derived.  The  derivation  of  the  stability  condition  by 
the  probabilistic  method  for  the  general  FDMTH  for 
transient  heat  conduction (diffusion)  is  included  as 
Appendix  C.  The  stability  condition  derived  by  the 
probabilistic  method  is  the  same  as  the  stability 
condition  derived  by  the  matrix  and  coefficient  methods. 

Summary. 

All  methods  of  stability  analysis  are  equivalent 
in  that  each  lead  to  the  same  stability  condition. 


< 


aAt  [ 


(1"fx)  <1“fy> 

- -  +  _ ] 

(Ax)2  (Ay)2 


<  h 


(33) 
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V.  Procedures 


Computer  System. 

The  data  in  this  study  is  from  programs  executed 
on  the  Harris  500  computer  of  the  Air  Force  Institute 
of  Technology. 

Computer  Programs. 

Programs  EXPLI,  IMPLI,  and  PRADI  were  written  for 
execution  of  the  fully  explicit ,  fully  implicit  and  Crank 
Nicholson  implicit,  and  Peaceman-Rachford  ADI  FDMTHs, 
respectively.  The  Guass-Sidel  method  was  used  in  IMPLI 
and  the  Thomas  method  was  used  in  PRADI.  Band  storage, 
see  Appendix  F,  was  used  to  minimize  memory  requirements. 

Error  Analysis. 

The  four  FDMTHs  studied  were  compared  on  the  basis 
of  truncation  error  and  stability.  Truncation  error  was 
computed  as  the  difference  between  the  FDMTH  temperature 
and  the  analytic  temperature.  The  analytic  temperature 
was  computed  from  the  analytic  solution  derived  at 
Appendix  A.  Stability  was  studied  by  counting  the  number 
of  oscillations  about  the  true  solution  and  noting  occur¬ 
rences  of  instability.  Only  data  for  those  nodes  along 
the  diagonal  from  the  northwest  corner  of  the  rectangular 
region  to  its  center,  see  Figure  13,  were  considered. 
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(*  -  nodes  for  which  data  is  collected) 


Figure  13.  Nodes  for  which  Data  Would  be  Collected 
for  a  11  by  11  nodal  array. 


Stability  Conditions. 

The  stability  condition  for  the  general  FDMTH  for 
two-dimensional  transient  heat  conduction (diffusion)  was 
derived  by  the  matrix  method,  the  coefficient  method, 
and  the  probabilistic  method.  The  stability  condition 
was  used  to  select  the  program  input  parameters  for 
execution  of  the  various  FDMTHs . 


VI.  Results 


Most  of  the  data  collected  in  this  study  is  at 
Appendix  E.  The  following  summarizes  the  results  of  this 
study. 

Stability. 

Applicability  of  the  general  stability  condition, 
described  in  Section  IV,  was  verified  for  the  fully  explicit 
and  fully  implicit  FDMTHs .  Though  the  Crank -Nicholson 
implicit  and  Peaceman-Rachford  ADI  FDMTHs  were  expected 
to  be  unconditionally  stable,  some  unstable  oscillatory 
behavior  was  observed  for  large  time  steps.  Unstable 
fP  oscillations  for  the  fully  explicit  FDMTH  were  typically 

as  depicted  in  Figure  14.  Unstable  behavior  of  the 
Crank-Nicholson  implicit  and  Peaceman-Rachford  ADI  FDMTHs, 
for  large  time  steps,  is  illustrated  by  figures  15  and 
16  through  18,  respectively.  Edge  effects,  stable  oscil¬ 
lations  for  nodes  near  corners,  were  observed  for  all 
FDMTHs.  Edge  effects  for  the  fully  implicit  FDMTH  are 
represented  by  Figure  19,  where  the  number  of  oscillations 
for  each  node  is  represented  by  a  number  in  the  same 
relative  position  as  the  position  of  the  node  in  the 
nodal  array. 
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X,  DISTANCE  FROM  THE  WEST  FACE  OF  THE  RECTANGULAR  REGION 

Figure  14.  Unstable  Oscillations  for  the  Fully  Explicit 
FDMTH  for  a  21  by  21  Nodal  Array.  Elapsed 
Time  is  0.08.  Time  Step  is  0.01.  Number  of 
Iterations  is  8.  tc  =  8. 


TEMPERATURE  (10 


Analytic  Solution 


r  , 


X,  DISTANCE  FROM  THE  WEST  FACE  OF  THE  RECTANGULAR  REGION 

Figure  15.  Unstable  Oscillations  for  the  Crank- 

Nicholson  Implicit  FDMTH  for  a  11  by  11 
Nodal  Array.  Elapsed  Time  is  0.3.  Time 
Step  is  0.1.  Number  of  Iterations  is  3. 
ic  =  10. 
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Discretization  Error . 


The  Crank-Nicholson  implicit  FDMTH  has  the  smallest 
discretization  error  for  large  numbers  of  time  steps. 

The  fully  explicit  FDMTH  is  the  next  most  accurate  FDMTH 
for  large  numbers  of  time  steps.  The  fully  implicit 
FDMTH  has  a  positive  discretization  error  for  large 
numbers  of  time  steps  that  is  larger  than  the  discret¬ 
ization  error  for  the  Crank-Nicholson  implicit  FDMTH. 

For  a  sufficiently  small  time  step,  the  discretization 
errors  of  the  fully  explicit  and  fully  implicit  FDMTH s 
are  nearly  equal.  For  large  numbers  of  time  steps, 
the  Peaceman-Rachford  ADI  FDMTH  is  the  least  accurate 
of  the  FDMTH s  studied.  The  Peaceman-Rachford  ADI  FDMTH 
has  a  negative  discretization  error  that  not  significantly 
improved  by  either  decrease  in  time  step  size  or  increase 
in  nodal  density.  Figure  20  is  a  comparison  of  all  FDMTH s 
studied  for  a  time  step  size  of  0.0025  and  an  elapsed 


time  of  1. 
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Execution  Time. 


Data  on  the  execution  times  of  the  different  FDMTHs 
was  not  collected.  In  general,  the  fully  implicit  and 
Crank -Nichols on  implicit  FDMTHs  were  slower  because  of 
the  use  of  the  Guass-Sidel  method  of  solution.  Figures 
21  and  22  depict  the  number  of  Guass-Sidel  iterations 
required  for  the  fully  implicit  and  Crank-Nicholson 
implicit  FDMTHs,  respectively,  for  a  11  by  11  nodal  array 
and  differing  size  and  number  of  time  steps. 


ELAPSED  TIME 


Number  of  Guass-Sidel  Iterations  Required  per 
Time  Step  for  the  Fully  Implicit  FDMTH  for  a 
11  by  11  Nodal  Array. 


Number  of  Guass-Sidel  Iterations  Required  per 
Time  Step  for  the  Crank-Nicholson  Implicit 

FDMTH  for  a  11  by  11  Nodal  Array. 


VII.  Discussion 

All  FDMTHs  considered  gave  unrealistic  results  for 
a  time  step  of  1  when  the  maximum  initial  temperature# 
TINIT  in  programs  EXPLI,  IMPLI,  and  PRADI;  a,  the  thermal 
diffusivity;  and  the  dimensions  of  the  rectangular  region 
were  equal  to  1.  The  order  of  the  truncation  error  for 
all  FDMTHs,  per  Table  2,  is  also  1  when  the  size  of  the 
time  step  is  1.  The  observed  unrealistic  results  for 
a  time  step  of  1  are  probably  due  to  discretization  error 
rather  than  the  type  of  numerical  instability  addressed 
by  the  stability  condition. 

The  accuracy  of  all  FDMTHs  improved  with  decrease  in 
the  size  of  the  time  step  and/or  increase  in  nodal  density 
The  accuracy  of  FDMTHs,  relative  to  each  other,  remained 
as  suggested  by  the  order  of  the  truncation  error  as  long 
as  the  size  of  the  time  step  was  small  enough  for  the 
FDMTHs  considered  to  remain  stable.  The  Peaceman-Rachford 
ADI  FDMTH  was  the  least  accurate  of  the  four  FDMTHs 
considered  for  large  numbers  of  time  steps  when  compared 
to  stable  FDMTHs  using  a  common  size  of  time  step  and 
nodal  array  size. 

Figures  2l  and  22  suggest  that  the  fully  implicit 
and  Crank-Nicholson  implicit  FDMTHs  require  fewer  Guass- 
Sidel  iterations  per  time  step  for  smaller  time  step  sizes 


VIII.  Conclusions 


Stability. 

The  general  stability  condition  correctly  predicts 
the  onset  of  instability  for  the  fully  explicit  FDMTH. 

The  fully  implicit  FDMTH  is  always  stable,  as  expected. 
The  Crank-Nicholson  implicit  and  Peaceman-Rachford  ADI 
FDMTHs,  expected  to  be  unconditionally  stable,  are 
unstable  for  large  time  step  sizes  for  the  initial¬ 
boundary  value  problem  studied. 

Accuracy. 

For  sizes  of  time  step  satisfying  the  general 
stability  condition  for  all  four  of  the  FDMTHs  considered 
for  large  numbers  of  time  steps,  the  Crank-Nicholson 
implicit  FDMTH  is  the  most  accurate  FDMTH  and  the 
Peaceman-Rachford  ADI  FDMTH  the  least  accurate  for  the 
initial-boundary  value  problem  studied. 

Complexity. 

The  simplest  FDMTH  algorithm  was  that  of  the  fully 
explicit  FDMTH. 

Stability  Condition. 

Derivation  of  the  general  stability  condition  by 
either  the  matrix  method,  coefficient  method,  or  proba- 


bilistic  method  leads  to  the  same  result, 

U-f  )  (1-f  ) 

■c  =  aAt [ - ]  <  H  (34) 

(Axr  (Ay) z 

Further,  while  a  FDMTH  may  be  stable  when  the  general 
stability  condition  predicts  instability,  the  accuracy 
of  a  FDMTH,  such  as  the  Peaceman-Rachford  ADI  FDMTH,  is 
significantly  improved  if  the  size  of  the  time  step  is 
one  for  which  the  general  stability  condition  predicts 
stability. 


IX.  Recommendations 


Use  of  This  Thesis. 

This  thesis  should  be  used  as  a  reference  or  starting 
point  for  future  studies  of  FDMTHs  to  encourage  the  study 
of  three-dimensional  FDMTHs  and  FDMTH  modelling  of  nuclear 
effects. 

Follow-on  Studies. 

A  follow-on  study  is  needed  to  better  define  the 
regions  of  stability  for  two-dimensional  FDMTHs,  possibly 
in  a  manner  similar  to  Table  3  and  Figure  12  of  this  study. 

Another  study  should  consider  FDMTH  solution  of  a 
different  initial-boundary  value  problem  having  a  damped 
analytic  solution. 

The  use  of  graded  spatial  grids  should  be  considered 
as  a  method  of  reducing  edge  effects. 
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APPENDIX  A 


Analytic  Solution  of  the  Initial-Boundary  Value  Problem  for 
Two-Dimensional  Transient  Heat  Conduction (Diffusion) 

An  initial-boundary  value  problem  describing  two- 
dimensional  transient  heat  conduction (diffusion)  for  a 
rectangular  region  is 

3T  3  2T  32T 

—  =  a  [  -  +  -  ]  (A-la) 

3t  3x2  3y2 

T(0,y,t)  =  T (a,y,t)  =  T(x,0,t)  =  T(x,b,t)  =  0  (A-lb) 

T(x,y,0)  =  simrxsimry  (A-lc) 

where  T  is  for  temperature,  a  is  the  thermal  diffusivity, 
x  and  y  are  spatial  dimensions,  and  t  is  for  time. 

The  solution  for  T(x,y)  follows. 

Let  T (x,y,t)  =  F(x,y)G(t) , 

then  3T/3t  =  T.  =  FG. ,  32T/3x2  =  T  -  F  G, 

L  U  Aa  AA 

and  32T/32y  =  T^  =  FyyG. 

Now,  FGt  =  a[FxxG  +  FyyG] , 


or  G./G  ■  X  ■  a[F  /F  +  F„/F] 


Let  F(x,y)  =  H(x)Q(y),  then  Fxx  =  HxxQ,  and  Fyy  =  HQ^  . 

x  =  a[HxxQ/HQ  +  HQyy/HQ]  =  a [Hxx/H  +  Qyy/Q] 

X/a  =  H  /H  +  Q  /Q 
xx  yy 

H  /H  ■  £  ■  X/a  -  Q  /Q 
xx  yy 

Solve  H  /H  =  £  for  H. 
xx 

[D2  -  £]H  =  0  H  -  fiisinh/^fx*  =  ezCOsh/^^T. 

Applying  the  boundary  condition  H(0)  *  0  , 

0  =  &2cosh0  -*■  02  =  0,  or  H  =  6isinh/-£x‘ , 

or  H  *  cisin/fx1  (by  the  trigonometric  identity: 

sinh/-l  u  —  sin/u7). 

4 

Applying  the  boundary  condition  H(a)  -  0  , 

0  =  fiisin/fa’  -*■  imr  «  a/f’  K  =  nt2ir2/a2,  m=0,l,2,  ... 
or 

H  =  Cisin[mirx/a]  . 

Similarly,  Q  =  disin[nny/b]  , 

or  F  *  Ci&isin[mirx/a] sin[mry/b]  . 
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When  solving  for  Q,  it  is  also  found  that 

nir  =  /X/a  -  |b  -*•  n2ir2/b2=  X/a  -  £  , 

or  X  =  a{n2u2/b2  +  £)  =  a(n2ir2/b2  +  m2TT2/a2)  . 

Solving  Gt/G  =  X,  for  G, 

[D  -  X]G  =  0  -*■  G  =  fexp[a  (m2ir2/a2  +  n2,rr2/b2)]t 

or, 

T(x,y,t)  =  C  exp[  -a  (m2ir2/a2  +  n2ir2/b2)  t]  sin[mirx/a]  sintniTy/b] 
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APPENDIX  B 


Matrix  Method  Derivation  of  the  Stability  Condition  for 
the  Two-Dimensional  Finite  Difference  Approximation  of 
Transient  Heat  Conduction (Diffusion) 


The  general  two-dimensional  finite  difference 
approximation  of  transient  heat  conduction (diffusion)  is 


Tk+1  _ 

P  P 


cxAt 


(Ax) 


— [  f  T^+^  +  (1-f  )T^ 
2 1  EE  1  rEl  E 


+  fWTW  +  (1_fw)Tw 


-  2f  Tp+1  -  2(l-f_  )T*] 

EW  r  EW  r 


aAt 


(Ay) 


fA+1  *  <1-fN)TN 


+  f_T*+1  +  ( 1-f c) T^ 


or  as  depicted  on  the  next  page. 


T(0,y)  =  T(l,y)  =  T  (x,  0)  =  T  (: 


Solving  for  T, 


T(x,y)  =  F (x) G (y) 


T  =  F  G 
XX  XX 


T  =  FG 

yy  yy 


or  X  FG  =  F  G  +  FG 

s  xx  yy 


X  -  F 

S  5 


-F 


XX 


+  X 


F 


c  = 

G 


-F  +  (X  -£)F  =  0 
xx  s 

[D2  +  (£-X  )]F  =  0  ,  let  n  = 
s 

F  =  eiSinhZ-nx'  +  fi2  sinh /-nx1 

F  ( 0 )  =  0  =  62cosh0  •>  62  =  0 

F  =  61  sinh/-rix‘ 

F  (1)  =  0  =  Cisin/rT  -*■  /rT  =  st 
F  =  Cisin(sirx) 


for  Gyy  -  £G  =  0,  as  for  F , 


G  =  aisin(tiry) 


and  T  «  C  sin(sTrx)  sin(sny) 


Let  T  =  sin(sirx)  sin(tiry)  . 


A  finite  difference  representation  of  equation  (B-l) 


b'T(i+l, j)  -  2b'T(i,j)  +  b'T(i-l, j) 

+  c'T(i, j+1)  -  2c'T(i,j)  +c'T(i,j-l)  =  X  T(i,j)  (B-6) 


or  as  shown  on  the  next  page. 
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mirsi  mrrs  mrtj 

b'sin[  -  +  -  3  sin[  -  ] 

MM  N 


rtMTsi 

-  2b1 sin [  - 

M 


nut  j 

sin[  - 

N 


mrrsi  mirs 

+  b'sin[  -  -  -  ]  sin[ 

M  M 


mirsi  mrtj  nirt 

c'sin[  -  j  sin[  -  +  -  ] 

M  N  N 


irnrsi  mrtj 

2c' sin [  -  ]  sin[  - 

M  N 


nnrsi 

+  c'sin[  - 

M 


nirtj 

sin[  - 

N 


irnrsi  mrtj 


X  sin[ 


1  sin  [ 


3 


Using  the  trigonometric  identity, 

sinacosB  =  *ssin(a+B)  +  *ssin(a-B)  , 

and  dividing  through  by  sin (musi/M) sin (mrt j/N) ,  equation 
(B-7)  reduces  to 

2b'cos  [imts/M]  -  2b'  +  2c 1  cos  [n-rrt/N]  -  2c'  =  X  (B-8) 

s 

By  the  trigonometric  identity, 

2 

sin  a  =  %(l-cos2a)  , 
equation  (B-8)  reduces  to 

X  *=  -4  [  b' sin2  (imts/2M)  +  c '  sin2  (nnt/2N)  (B-9) 

s 

For  equation  (B— 9 )  to  be  stable,  must  have  that 

1  <_  |  1  -  4[  b' sin2  (mus/2M)  +  c' sin2  (nirt/2N)  ]  | 

(B-10) 

or 

-1  <  1  -  4[  b' sin2  (mirs/2M)  +  c' sin2 (mit/2N)  ]  (B-ll) 
.  2 

or,  since  sin  a  <  1, 

h  >  b'  +  c’  .  (B-12) 

Or,  the  general  stability  condition  is 

(1-f  )  (1-f  ) 

tc  =  aAt  l  - ^  ^  ]  <  h  (B-13) 

(Ax)  z  (Ay)  £ 
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for  the  Two-Dimensional  Finite  Difference  Approximation  of 


Transient  Heat  Conduction (Diffusion) 


The  general  two-dimensional  finite  difference 
approximation  of  transient  heat  conduction (diffusion)  is 
as  given  in  Appendix  B,  equation  (B-l) ,  or 

( l+2p+2q) Tp+1  -  PT*+1  -  Pt£+1  -  qT*+1  -  qT*+1 

-  (l-2p,-2ql )Tp  +  p 'Tg  +  P't£  +  q'T*  +  q'T^  (C- 


where  p  =  aAtfx/ (Ax) 
q  =  aAtfy/(Ay)2 
p'  «  aAt(l-fx)/(Ax) 
q'  =  aAt (1-fy) / (Ay) 

fx  =  fE  =  fW  =  fPmj 


and  f  =  f  =  f  =  f 

y  N  b  *NS 

k+1  k 

The  ratio  of  Tp  to  Tp  ,  X,  from  equation  (C— 1 )  is 

(l-2p'-2q')  p(T*+1  +  T*+1)  +  q(T*+1  +  T**1) 


(l+2p+2q) 


(l+2p+2q) 


,  P’ (IE  +  TW>  *  q' (TS  +  TSI  (c.2 

(l+2p+2q) 


Both  Myers  and  Patankar  consider  only  the  case  of  a  single 
interior  node  so  that  T£  =  Tw  =  TN  =  Tg  =  0.  For  a  single 
interior  node. 


X  = 


l-2p'-2q' 

l+2p+2q 


(C-3) 


The  stability  curves  predicted  by  (C-3)  are  depicted 
in  Figure  C-l.  The  general  stability  condition  can  be 
derived  from  equation  (C-3)  as  follows: 


1  -  2<xAt  (l-f„)  /  (Ax)  2  -  2aAt (1-f  )  /  (Ay)  2 


1  -  2aAtf  /(Ax)2  +  2aAtf  /(Ay)2 
x  y 


(C-4) 


For  stability,  X  must  be  positive.  The  denominator  of 
equation  (C-6)  is  always  positive,  so  only  the  numerator 
need  be  considered  further.  For  stability  then 

0  <  1  -  2aAt (1-f  ) / (Ax) 2  -  2aAt(l-f  )/(Ay)2 

*  x  y 


h  >  aAt[ (1-f  ) / (Ax) 2  +  (1-f  ) / (Ay) 21  . 
x  y 


Or,  the  general  stability  condition  is 


k  =  aAt [ 


(1-f  )  (l-f„) 


1  <  %  • 


(C-5 


Probabilistic  Method  Derivation  of  the  Stability  Condition 
for  the  Two-Dimensional  Finite  Difference  Approximation  of 

Transient  Heat  Conduction (Diffusion) 

The  general  two-dimensional  finite  difference  approxi¬ 
mation  of  transient  heat  conduction (diffusion)  is  as  given 

k+1 

in  Appendix  B,  equation  (B-l) .  Let  Pe  be  the  probability 
of  node  P  attaining  the  potential (temperature)  of  node  e 
at  time  t  ±  iAt,  then 


,k+l 

P 

=  1 

(D-l) 

,k+l 

E 

=  aAtfE/ (W) (Ax) 2 

(D-2) 

,k+l 

W 

=  aAtfw/(W) (Ax)2 

(D-3) 

,k+l 

N 

=  aAtfN/(W) (Ay)2 

(D-4) 

,k+l 

S 

=  aAtfg/ (W) (Ay) 2 

(D-5) 

0 

0- 

[1  -  2aAt (l-f_  ) / (Ax) 2  -  2aAt(l-f„  )/(Ay)2] 


(D-ll) 


and  W  =  1  +  aAtf  /(Ax)  +  aAtf  / (Ay) *  . 

EW  PNS 

Using  the  probabilistic  coefficients  defined  on  the 

previous  page,  the  general  FDMTH  becomes 

k+1  k+1  _  k+1  k+1  pk+l_k+l  k+1  k+1  k+1  k+1 

PP  TP  “  PE  TE  +  PW  TW  +  PN  TN  +  PS  TS 

+  PETE  +  PWTW  +  PK  +  PSTS  +  PPTP  •  (D‘12> 

As  long  as  all  f  are  less  than  or  equal  to  1,  the  only 

probabilistic  coefficient  that  might  be  negative  is  Pp  . 
Since  negative  probabilities  are  not  realistic,  the 
following  stability  condition  is  inferred. 


1  -  2a At [ 


(i-fj  u-f  ) 


(Ax) 


(Ay) 


1  >  0 


(D-13) 


or 


_  ...  (1"fx>  . 

k  =  aAt  [  - =■  +  *  ^ 


(Ax) 


(Ay) 


1  <  h  , 


(D-14' 


where 


E 


EW 


f  c  f 

y  n 


'NS 


Equation  (D-14)  is  the  general  stability  condition  for 
the  two-dimensional  finite  difference  approximation  of 
transient  heat  conduction (diffusion) . 


APPENDIX  E 


FDMTH  Data 

The  volume  of  data  generated  in  this  study  was  too 
great  to  graphically  protray  all  of  it.  The  data  presented 
in  this  appendix  represents  most  of  the  more  significant 
data  collected  and  not  previously  described  in  Section  VI, 
Results . 


E-l-1  Unstable  Oscillations  for  the  Peaceman- 
Rachford  ADI -FDMTH  for  a  31  by  31  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  10.  Time  Step 
is  1.  Number  of  Iterations  is  10. 

E-l-2  Unstable  Oscillations  for  the  Peaceman- 
Rachford  ADI  FDMTH  for  a  31  by  31  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  20.  Time  Step 
is  1.  Number  of  Iterations  is  20. 

E-2-1  Comparison  of  FDMTH s  for  a  11  by  11  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  6.  Time  Step  is 
1.  Number  of  Iterations  is  6. 

E-2-2  Comparison  of  FDMTHs  for  a  11  by  11  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  0.2.  Time  Step 
is  0.01.  Number  of  Iterations  is  20. 

Comparison  of  FDMTHs  for  a  11  by  11  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  1.  Time  Step  is 
0.001.  Number  of  Iterations  is  1000. 


E-2-3 


Figure 


E-2-4 


E-2-5 


Page 


E-2-6 


E-3-1 


E-3-2 


Comparison  of  FDMTHs  for  a  11  by  11  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  0.1.  Time  Step 
is  0.0001.  Number  of  Iterations  is  1000. 

Comparison  of  FDMTHs  for  a  21  by  21  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  2.  Time  Step  is 
1.  Number  of  Iterations  is  2. 

Comparison  of  FDMTHs  for  a  21  by  21  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  0.01.  Time  Step 
is  0.0001.  Number  of  Iterations  is  100. 

Number  of  Oscillations  by  Node  for  the 
Peaceman-Rachford  ADI  FDMTH  for  a  11  by  11 
Nodal  Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  0.5.  Time  Step 
is  0.0001.  Number  of  Iterations  is  5000. 

Number  of  Oscillations  by  Node  for  the 
Crank-Nicholson  Implicit  FDMTH  for  a  21  by 
21  Nodal  Array  Imposed  Over  a  1  by  1 
Rectangular  Region.  Elapsed  Time  is  0.01. 
Time  Step  is  0.0001.  Number  of  Iterations 
is  100. 


The  following  acronyms  are  used  in  the  figures  of  this 
appendix . 

ANALY  Analytic  Solution 

EXPLI  Fully  Explicit  FDMTH 

IMPLI  Fully  Implicit  FDMTH 

CNICH  Crank-Nicholson  Implicit  FDMTH 


PRADI 


Peaceman-Rachford  ADI  FDMTH 


TEMPERATURE 


Analytic  Solution 
is  approximately  zero 


X,  DISTANCE  FROM  THE  WEST  FACE  OF  THE  RECTANGULAR  REGION 


Figure  E-l-2 


Unstable  Oscillations  for  the  Peaceman- 
Rachford  ADI  FDMTH  for  a  31  by  31  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  20.  Time  Step 
is  1.  Number  of  Iterations  is  20. 

900. 


71 


TEMPERATURE 


TEMPERATURE (10 


TEMPERATURE (10 


I 


X,  DISTANCE  FROM  THE  WEST  FACE  OF  THE  RECTANGULAR  REGION 


Figure  E-2-3.  Comparison  of  FDMTHs  for  a  11  by  11  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  1.  Time  Step  is 
0.001.  Number  of  Iterations  is  1000. 
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TEMPERATURE 


i(il 


LEGEND 

_ 

ANALY 

EXPLI 

IMPLI 

CNICH 

PRADI 

X,  DISTANCE  FROM  THE  WEST  FACE  OF  THE  RECTANGULAR  REGION 

Figure  E-2-4.  Comparison  of  FDMTHs  for  a  11  by  11  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  0.1.  Time  Step  is 
0.0001.  Number  of  Iterations  is  1000. 


TEMPERATURE 


LEGEND 


-O-  ANALY 
~0~  EXPLI 
-Or  CNICH 
PRADI 
IMPLI 


X,  DISTANCE  FROM  THE  WEST  FACE  OF  THE  RECTANGULAR  REGION 


Figure  E-2-6 


Comparison  of  FDMTHs  for  a  21  by  21  Nodal 
Array  Imposed  Over  a  1  by  1  Rectangular 
Region.  Elapsed  Time  is  0.01.  Time  Step 
is  0.0001.  Number  of  Iterations  is  100. 
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is  5000 


Step  is  0.0001.  Number  of  Iterations  is  100 


APPENDIX  F 


Use  of  Band  Storage  to  Reduce  Computer  Program 
Memory  Requirements  for  Coefficient  Matrices* 

This  appendix  presents,  by  example,  the  use  made  of 
band  storage  for  storage  of  banded  coefficient  matrices 
characteristic  of  the  fully  implicit,  Crank-Nicholson 
implicit,  and  Peaceman-Rachford  alternating  direction 
implicit  (ADI)  finite  difference  methods  (FDMTHs) . 

For  an  m  by  n  nodal  array  imposed  over  a  rectangular  region, 
application  of  any  of  the  forementioned  FDMTHs  results  in 
a  banded  coefficient  matrix.  For  the  Peaceman-Rachford  ADI 
FDMTH,  the  resulting  coefficient  matrix  is  tridiagonal. 

For  the  fully  implicit  and  Crank-Nicholson  implicit  FDMTHs, 
the  resulting  coefficient  matrix  has  a  band  width  of  2M-1 
or  2N-1  depending  on  the  order  in  which  the  nodes  are 
numbered.  M  =  m  -  1  and  N  *  n  -  1.  The  variables  M  and 

N  are  introduced  for  a  convenience  in  coefficient  subscript- 

✓ 

ing  described  later.  The  remainder  of  this  appendix  is 
concerned  with  the  coefficient  matrices  of  the  fully 
implicit  and  Crank-Nicholson  implicit  FDMTHs.  The  coeffi¬ 
cient  matrix  of  the  Peaceman-Rachford  ADI  FDMTH  may  be  band 
stored  with  minor  modifications  to  the  following. 


*Alan  Jennings,  Matrix  Computations  for  Engineers  and 

Scientists ,  (London:  John  Wiley  and  Sons,  Inc.,  1977) 

pp.  95-96. 
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The  coefficient  matrix  is  dimensioned  MDL  by  MDL 
for  the  fully  implicit  and  Crank-Nicholson  FDMTHs  where 
MDL  =  (M-l)  times  (N-l) .  In  the  nodal  array  is,  for 
example,  31  by  31,  then  M  =  N  =  30  and  the  coefficient 
matrix  is  841  by  841.  By  taking  advantage  of  the  banded 
nature  of  coefficient  matrices  for  the  fully  implicit  and 
the  Crank-Nicholson  FDMTHs;  and,  noting  that,  when  using 
the  Guass-Sidel  solution  technique,  no  use  is  made  of  the 
zero-valued  coefficients  not  within  the  banded  portion  of 
the  coefficient  matrix;  the  coefficient  matrix  can  be 
represented  by  a  band  storage  array  dimensioned  (2M-1)  by 
MDL.  For  a  31  by  31  nodal  array,  the  band  storage  array 
is  dimensioned  59  by  841,  or  less  than  one  tenth  of  the 
size  of  the  coefficient  matrix  that  it  replaces,  dimensioned 
841  by  841.  The  following  example  illustrates  the 
application  of  band  storage. 

Consider  the  coefficient  matrix  for  the  fully  implicit 
FDMTH  for  a  5  by  5  nodal  array,  shown  in  Figure  F-l,  where 
coefficients  not  shown  are  zero-valued  and  outside  the 
band  of  interest.  In  this  case,  M  =  N  =  4  and  the  band 
width  is  7.  If  zeros  are  added  to  the  band  of  Figure  F-l 
so  that  the  band  resembles  a  parallelogram  as  shown  in 
Figure  F-2,  then  a  new  coefficient  matrix  of  dimension 
7  by  9  can  be  defined  whose  nonzero  elements  would  be  as 
depicted  in  Figure  F-3. 
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Figure  F-l.  A  Coefficient  Matrix  for  a  5  by  5  Nodal  Array 


Figure  F-2.  Augmented  Coefficient  Matrix  for  a  5  by  5 
Nodal  Array. 
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Figure  F-3.  Coefficients  of  Interest  in  Band  Storage 


The  form  of  band  storage  represented  by  Figure  F-3 
is  perferred  because  it  maintains  the  node  to  matrix  row 
subscript  relationship  (i.e.  coefficients  b^  result  from 
application  of  a  FDMTH  at  node  4) .  A  general  form  of  the 
band  stored  coefficient  matrix  is  as  depicted  in  Figure 
E-4. 

Figure  F-4  also  illustrates  the  reason  for  introducing 

the  variables  M  and  N.  Without  M  a'.d  N,  the  coefficient 

b-  „  would  become  b~  „  .  and  coefficients  on  the  main 
1,M  2,m-l 

diagonal  would  become  bm_^  ^ ,  where  i  =  1,2,  ...  .  M  and 
N  are  used  because  the  resulting  subscripts  of  the  band 
are  preferred. 


1  ,MDL  M-2 ,MDL  M- 1 , MDL  M,MDL 
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